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Abstract. Motivated by applications to numerical simulation of flows in highly het- 
erogeneous porous media, we develop multiscale finite element methods for second 
order elliptic equations. We discuss a multiscale model reduction technique in the 
framework of the discontinuous Galerkin finite element method. We propose three dif- 
ferent finite element spaces on the coarse mesh. The first space is based on a local 
eigenvalue problem that uses a weighted L2— norm for computing the "mass" matrix. 
The second space is generated by amending the eigenvalue problem of the first case 
with a term related to the penalty. The third choice is based on generation of a large 
space of snapshots and subsequent selection of a subspace of a reduced dimension. 
The approximation with these spaces is based on the discontinuous Galerkin finite el- 
ement method framework. We investigate the stability and derive error estimates for 
the methods and further experimentally study their performance on a representative 
number of numerical examples. 

1. Introduction 

In this paper we present a study of numerical methods for the simulation of flows 
in highly heterogeneous porous media. The media properties are assumed to contain 
multiple scales and high contrast. In this case, solving the systems arising in the ap- 
proximation of the flow equation on a fine-grid that resolves all scales by the finite 
element, finite volume, or mixed FEM could be prohibitively expensive, unless special 
care is taken for solving the resulting system. A number of techniques have been pro- 
posed to efficiently solve the these fine-grid systems. Among these are multigrid methods 
(e.g., OE]), multilevel methods (e.g., [SI1E2]), and domain decomposition techniques 

(e.g., [m [IHl [231 12a 1251 EO]). 

More recently, a new large class of accurate reduced-order methods has been intro- 
duced and used in various applications. These include Galerkin multiscale finite elements 
(e.g., la m tta [201 1211 122]), mixed multiscale finite element methods (e.g., P [21 [I [2Z] ) , 
the multiscale finite volume method (see, e.g., [2H]), mortar multiscale methods (see 
e.g., ini [33]), and variational multiscale methods (see e.g., [26]). Our main goal is to 
extend these concepts and develop a systematic methodology for solving complex multi- 
scale problems with high-contrast and no-scale separation by using discontinuous basis 
functions. 
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In this paper, we study the multiscale model reduction techniques within discontinuous 
Galerkin framework. As the problem is expected to be solved for many input parameters 
such as source terms, boundary conditions, and spatial heterogeneities, we divide the 
computation into two stages (following known formalism [HI [29]): offline and online, 
where our goal in the offline stage is to construct a reduced dimensional multiscale space 
to be used for rapid computations in the online stage. In the offline stage [I9], we generate 
a snapshot space and propose a local spectral problem that allows selecting dominant 
modes in the space of snapshots. In the online stage use the basis functions computed 
offline to solve the problem for current realization of the parameters (a further spectral 
selection may be done in the online step in each coarse block). As a result, the basis 
functions generated by coarse block computations are discontinuous along the coarse- 
grid inter-element faces/edges. Previously, e.g. [19], in order to generate conforming 
basis functions, partition of unity functions have been used. However, this procedure 
modifies original spectral basis functions and is found to be difficult to apply for more 
complex fiow problems. In this paper, we propose and explore the use of local model 
reduction techniques within the framework of the discontinuous Galerkin finite element 
methods. 

We introduce a Symmetric Interior Penalty Discontinuous Galerkin (SIPG) method 
that uses spectral basis functions that are constructed in special way in order to reduce 
the degrees of freedom of the local (coarse-grid) approximation spaces. Also we discuss 
the use of penalty parameter in the SIPG method and derive a stability result for a 
penalty that scales as the inverse of the fine-scale mesh. We show that the stability 
constant is independent of the contrast. The latter is important as the problems under 
consideration have high contrast. 

We also derive error estimates and discuss the convergence issues of the method. Ad- 
ditionally, the efficacy of the proposed methods is demonstrated on a set of numerical 
experiments with fiows in high-contrast media where the permeability field has subre- 
gions of high conductivity, which form channels and islands. In both cases we observe 
that as the dimension of the coarse-grid space increases, the error decreases and the 
decrease is proportional to the eigenvalue that the corresponding eigenvector is not in- 
cluded in the coarse space. In particular, we present results when the snapshot space 
consists of local solutions. 

The paper is organized in the following way. In Section |2] we present our model problem 
in a weak form and introduce the approximation method that involves two grids, fine 
(that resolves all scales of the heterogeneity) and coarse (where the solution will be 
sought). On each cell of the coarse mesh we introduce a lower dimensional space of 
functions that are defined on the fine mesh. We also show that the method is stable in a 
special DG norm. In Section |3] we present three different choices of local spaces. The first 
two are based on few eigenfunctions of special spectral problems in the style of [HI [23] . 
The third choice is based on the concept of snapshots, e.g. [8|l29]. In Section|l]we present 
some numerical experiments and report the error of the Discontinuous Galerkin method 
with the constructed coarse-grid spaces and in Section[5l we discuss the numerical results. 
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The theoretical results are derived under the assumption that the penalty stabilization 
depends on the fine-mesh size. Based on the numerical experiments we can conclude that 
the interior penalty Galerkin method gives reasonable practical results in using coarse- 
grid spaces generated by special problems, solved locally on each coarse-grid block, that 
take into account the highly heterogeneous behavior of the coefficient of the differential 
equation (in our case, the permeability). 

2. Continuous and discrete problems 
We consider the following problem: Find u* G Hq{Q) such that 
(1) a{u*,v) = f{v) hi aA\v e H^{n) 

where 

a{u,v) := Ik{x)Vu ■ Vvdx and f{v) := / fvdx. 
Jn Jq 

Here is a bounded domain in W^, d = 2,3 with polygonal boundary. We assume that 

/ G L2{^1) and the coefficient k,{x) represents the permeability of a highly heterogeneous 

porous media with high contrast, that is high ratio between the maximum and minimum 

values, see Figure [H Our main goal in this paper is to develop an approximation method 

for ([T]) on a coarse grid using certain "low energy" local eigenfunctions. 

We consider the two dimensional case. The method and results presented here extend 
for three dimensional case. We split the domain Q into disjoint polygonal subregions 
{^i}iLi of diameter 0{Hi) so that Q = U^ifli. We assume that the substructures 
{^i}iLi form a geometrically nonconforming partition of fl. In this case, for i ^ j, the 
intersection dQi fl dflj is either empty, a vertex of Qi and/or Qj, or a common edge of 
dQi and dQj. We recall that in the case of geometrically conforming decomposition, the 
intersection dQi (1 dflj is either empty or a common vertex of Qi and Qj, or a common 
edge of Qi and flj. Similar construction is assumed in 3-D with f2j being polyhedra. 

Further, in each Qi we introduce a shape regular triangulation Th{^i) with triangular 
elements and maximum mesh-size hi. The resulting triangulation of Q is in general 
nonmatching across dQi. Let Xh[Qi) be the regular finite element space of piecewise 
linear and continuous functions in Th{Qi). We do not assume that functions in Xh{^i) 
vanish on dQi fl dQ. We define 

and represent functions v of Xh{fl) as f = {vi}fLi with Vi G Xh{fti)- For simplicity, we 
also assume that the permeability k{x) is constant over each fine-grid element. 

Due to the fact that Th{^i) and Th{^j) are independent from each other on a common 
edge E = dfli fl dQj they may introduce two different partitions of E which are merged 
to obtain a set of faces Eij C E. Since the functions in Xh{^) are discontinuous along 
the interfaces, it is necessary to distinguish between E G Qi and E C Qj. From now 
on the fij-side of E will be denoted by E^j while on the f2j-side of E will be denoted by 
Eji. Geometrically, Eij and Eji are the same object. 
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We use the following harmonic averages along the edges Eij. For i,j G {1, . . . , A^} 
define 

Note, that the functions Kij and hij are piecewise constants over the edge E^j on a mesh 
that is obtained by merging the partitions Th{^i) and ThiS^j) along their common edge 

Eiy 

The discrete problem obtained by the DG method, see [3 [12] is: Find u\ = {u^ i}fLi ^ 
Xhi^l), ul^i e Xh{fli), such that 

(3) a^^K,t;^) = /K) for all vh = {vhAli ^ M^), 

where aj^^{uh,Vh), defined on Xh{^) x Xh{^), and /(f), defined on Xh{fi), are given by 

TV ^ /• 

(4) a^^(w, v) = Y, (^^''i^, ^) and f{v) = J2 f^idx. 

i=l i=l "^^i 

Here each local bilinear form af *^ is given as a sum of three symmetric bilinear forms: 

(5) af^iu, v) := ai{u, v) + Si{u, v) + pi{u, v), 
where Oj is the bilinear form associated with the "energy" , 

(6) ai{u,v) := / K{x)Vui ■ Vvidx, 
the Si is the bilinear form ensuring consistency and symmetry 

(7) Si{u, v):= J2 T I ""'^ ~ ""'^ + ^^""^ ~ 

and Pi is the penalty bilinear form that is added for stability 

1 (5 f 

(8) pi{u,v) := ^ TT~ I f^ijit^j -Ui){vj -Vi)ds. 

Here Kij is defined in ([2]) and denotes the outward normal derivative on dVli. The 
parameter 5 is a positive penalty parameter. In order to simplify notation we included 
the index j = 9 in the definition of the bilinear forms Si and Pi above. In order to include 
Eid '■= dQi n dfl in the summation sing, we set kj = 2 when i,j ^ d and kj = 1 when 
j = d. We also let vg = for all v G and define KiQ = and hiQ = hi. We note 

that when Kij is given by the harmonic average, then min{Kj, Kj} < Kij < 2min{Kj, Kj}. 
For later use we define the positive bilinear forms di as 

(9) di{u,v) = ai{u,v) + pi{u,v), 
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and the broken bilinear form dh for Xh{^)'- 

N 

(10) dhiu,v) := y^^di{u,v). 



i=l 



For u = {ui}]Li G Xh{^) the associated broken norm is then defined by 

N 



(11) II^IIm = MU,U) = ^l\\ K.!^Ui ||i2(f^^) + — — 1^ Kij{Ui-Ujfds ) . 

i=l I E.jCdQi 



We also have the following lemma shown in Lemma 3.1]. Here we provide a sketch 
of the proof for the sake of completeness. 

Lemma 2.1. There exists 5o > such that for 5 > 5q and for all u G X/i(fi) the 
following inequalities hold: 

Jodiiu, u) < af^{u, u) < 'jidi{u, n), i = l,...,N, 

and 

(12) -fodh{u,u) < a^^{u,u) < jidh{u,u), 

where 70 and 71 are positive constants independent of the Ki, hi Hi and u. 

Proof. First, we want to prove that '-yodi{u,u) < af^iu^u). Since af^{u,u) = ai{u,u) + 
Si{u,u) + di{u,u), the proof reduces to bound \si{u,u)\. Note that 



Si{u,u) 



El f 1 



EiiGdn, 



where we have defined Iij := Kij^{uj — Ui)ds. We have 



l/2v7 



L^(E„) 



l^ij [Uj - Ui) 



Using the following inequality for Ui G Xh{^t 

2 



hi 



V2v7 



L2(£;,,) 



< c 



l/2v7 



the Young's inequality with arbitrary e > and the fact hij < 2hi , we get 



< C 



k/ yui 



+ 



1 1 



1/2/ 



2 



Then, multiplying by l/kj and summing over the edges Eij C Qi, we get 
|si(M,M)| < 2CNEe 



V2v7 



I ^ sr^ 11 1/2/ N 



2CNEeai{u,u) + —pi{u,u). 

4eo 



C 



2 
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Here Ne denotes the number of edges of subdomain Vti. Choosing e = 1I{ACNe) we get 



and then 



C^Ne 

\si{u,u)\ < 0.5ai{u,u) H — Pi{u,u) 



C^Ne. , . , nn C^Nr 



0.5ai{u, u) + (1 — -^—r^)pi{u, u) < a^^{u, u) < 1.5aj(M, «) + (!+ ^ ^ )pi{u, u). 



Therefore the resuhs holds if we take 5 > 5o > CNe, 7o = min{0.5, 1 - {C'^Ne)/5} 
and 7i = max{1.5, {C'^Ne)/5}. □ 

Remark 2.2. We note that 71/70 in Lemma \2. 1\ deteriorates when 5 gets larger. In 
practice, however, 6 > 60 is chosen such that 6 = 0(1), therefore, from now on we 
assume that all the estimates will not depend on 6. 

3. Coarse-grid spaces 

In this section, we will construct local multiscale basis functions. We will follow 
GMsFEM where one needs the space of snapshots, see [171115]. In this space of snapshots, 
local spectral problems are designed and solved to compute multiscale basis functions. 
To keep our presentation simple, we first use the space of snapshots to be fine-grid 
functions within a coarse region. Thus, the local spectral problems will be posed on the 
fine grid. Next, we will discuss how a general space of snapshots can be used. 



3.1. Fine-grid snapshot space and weighted eigenvalue problem. Following [201 
[15] we consider the eigenvalue problem in Qi for the eigenvalues Af^ and the eigenf unc- 
tions ifjj^{x): 

(13) — div {k{x)\/'iI>1() = XliK,il>le, x G fij, k{x)ViI)1^^ ■ n = 0, x E dVLi, 

where n is the outer unit normal vector to dVti and k is a properly selected weight; for 
scalar permeability, we select k = k while for tensor permeability we refer to ^U\. The 
super-index / is used to distinguish from the other two methods we develop here (with 
indexes II and III). 

The eigenvalue problem considered above is solved in a discrete setting. Namely, for 
any given subdomain Qi find G Xh{Qi) such that 

(14) ai{iplf,, z) = Xlj,mi{iplf,, z) for all z G Xh{Qi). 
Here, aj is defined in ([6]) and the bilinear form mi{-,-) is defined by 



(15) 'mi{v, z) = Kvz. 

We order the eigenvalues so that < Af^ < Afg < •••• < ^1 n y where Ni is the number 
of vertices of Th{Qi), i.e., the number of degrees of freedom associated to Xh{Qi). Then, 
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in each subdomain we take the Li eigenfunctions corresponding to the smallest 
eigenvalues and use them as the multiscale basis. More precisely, define 

:= span 1 < / < L,} C X;,(a,), z = 1, . . . , iV. 

Finally, the coarse space is defined as 

Xj,{Q) := X^(fii) X ■ ■ ■ X X^(fijv) C Xhift). 

We refer to Xjj{fl) as an spectral coarse space due to its construction. 
Now the coarse-grid problem is: find m^*'^ G Xff{fl) such that 

(16) a^''(«r'^ vh) = f{vH) for all vh e Xh^. 

Note that the dimension of Xjj{Q) depends on the number of eigenvectors chosen in 
each coarse block Qj. An ideal situation would be when small number of eigenvectors in 
Qi represent (approximate) well the restriction of the solution to that subdomain. 

3.2. Fine-grid snapshot space with amended eigenvalue problem. Motivated 
by the error analysis developed below in Section 13.51 we use the following modified 
eigenvalue problem. Find ipj^ G Xh{^i) such that 

(17) a,(V'g, z) = \, {m,i^ll z) + mf (^-J^, z)) for all z G X^(a) 
where Oj is defined in (I6l), mj(-, ■) is defined by (ITSl) and 

1 (5 /* 

(18) mf{v,z)= TT~ I f^ijvzds. 

These eigenvalue problems allow us to obtain simple error estimates since the eigen- 
vectors can approximate fine functions simultaneously in a norm that includes interior 
weighted semi-norm in a coarse-grid block and weighted L^-norm on the interfaces. 

As before we order the eigenvalues as < Af^ < Af 2 < ■ ■ ■ < and we choose Li 
eigenfunctions corresponding to the smallest eigenvalues and use them as the multiscale 
basis. Define 

Xll{ni) := span {V^/j, 1 < / < C Xh{n^), i = l,...,N, 
and the coarse space 

x'Jin) := x'J{n,) X ■ ■ ■ X x'Jin^) c x^in). 

Now the coarse-grid problem similar to (fT6l) : find u^^'^^ G X'J (Q) such that 

(19) a^K^''', vh) = f{vH) for all vh G X^n). 
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3.3. General snapshot space and an example. In general, one can consider a gen- 
eral snapshot space for solving local eigenvalue problems. As we discussed in the In- 
troduction, the use of general snapshot space can have an advantage in case additional 
information is known about the local solution space. The subset of all possible function 
that satisfy the know properties of the unknown solution can be taken as the snapshot 
space. In this way we solve eigenvalue problem only on interesting (smaller dimension) 
subspaces instead of the space of fine degrees of freedom. For example, if solutions need 
to be computed only for a subspace of possible source terms, one can restrict the snapshot 
space to the space of local solutions for those sources and do not consider all fine-grid 
functions. To demonstrate that it is possible to use a snapshot space strictly smaller than 
Xh{Q), we consider an example where the snapshot space consists of all local solutions 
of the homogeneous equation with boundary conditions restriction on the boundary of 
the finite element nodal basis functions (or the set of all discrete Oj— harmonic functions 
in each block). More precisely, for the nodal basis function 6k{x) corresponding to the 
k—th node on dfli, we consider the problem 

(20) - div(KV0j,fc) = in Qi, (pi^k = h on dVLi, 

The (pi^k £ k = 1, . . . , Mi, is the (finite element) solution of this local problem. 

Here Mj denote the number of nodal basis function corresponding to nodes on dQi. Then 
the space of snapshots is defined by 

(21) ^P(a) = span{0,,fc, 1 < A; < M,}, t = l,...,N. 

Remark 3.1. Here, the reference solution we want to approximate on a coarse grid is 
the Galerkin projection ofu^, solution of into the global snapshot space 
Xf«^(fi,)x---xXr«^(fi^). 

Our objective is to construct a possibly smaller dimension space X^^^^(f2j) which is 
a subspace of X^^^^(fij). The construction is done using appropriate spectral decom- 
position. For this, define the matrices 

= [a.i<P^,k,M]k%=i, and Mf = [mf (0,,, 0,,O]5' 
and solve the following algebraic eigenvalue problem 

(22) = Aff PMf ^Pa,,. 

We write ai^i = (ai^i-i, . . . ,ai^i.Mi) ^ K*^' and define the corresponding finite element 
functions, ip-^ ^ ^h{^i) as 

Mi 

'^i/ = XI ^iAk'Pi,k, i=l,...,Mi. 

k=l 

Note that the matrices ^^^^P and Mf^^P are computed in the space of snapshots in 
Hi. Assume that < A^^^^P < ... < A^^^P, and choose the Li eigenvectors ^ • • • ) i^ij^i 
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that correspond to the smallest Li eigenvalues. We intoduce 

Xll\Vti) = span {V^g^ -.1 = 1,..., Li} for i = 1, . . . , iV 
and define the global coarse space by 

xjl\n) := x'J\n,) X ■ ■ ■ X xjl\n^) c x,(fi). 

The coarse problem is: find u^'^'^^^ G Xjj\Q) such that 

(23) «r vh) = f{vH) for all vh G Xll\n). 

Remark 3.2. Note that, according to the definition of mf in / figj) . the matrix m/'^^^ 
scales with 1 / hij . Then, the resulting eigenvalues scale with hij while the eigenspaces do 
not depend on hij . A similar situation is also valid for Method II and the eigenvalue 
problem [T7^ . It is easy to see from our main Theorem \3.5\ ( estated and proved below) that 
this scaling does not affect the convergence rate with respect to the number of eigenvectors 
used in the coarse space. 

Remark 3.3. Instead of the m/'^'^^ defined above, we can use 

M-^^P = [mi{(j)i^kAi,k') + mi{(f)i^k,(l)i,k')]^;k'=v 

3.4. Stability estimate. In this section, we present a best approximation result for 
the coarse-grid solution. 

Lemma 3.4. Let u\ G X^iyt) and u^^'^ G Xj^{Q) be the solutions of (QP and / f7^) . 
correspondingly. We have 

(24) dh{ul - u^''\ ul - u^'^) < Cidh{ul -v,ul- v) for all v G X^(^]) 
with Ci is independent of K,i,hi,Hi, and u^^'^ . 

Proof For all v G Xj^{n), 

(25) a^''iulv) = f{v) and a^«(«^^ t;) = /(t;). 
Then a^'^{u*j^ — uf^'^'\v) = and since u^^'^ G Xjj{Q), 

/■r,r\ DG f * msJ * ms,I\ DG f * msJ * \ 

Using a Cauchy-Schwarz inequality and f lT2]) in Lemma 12. 1^ 

7 / * ms,I * ms,I\ ^ DG / * msJ * ms,I\ DG t * msJ * \ 

^ DG I * msJ * msJ\\/2 DG f * * \ll2 

^ J / * ms,I * ms,I\l/2 J / * * \l/2 

Taking Ci = (71/70)^, we get 

7 / * ms.I * ms.I\ ^ J / * * \ 

and this completes the proof. □ 
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Analogous best approximation results, with respect to the fine-grid reference solutions, 
hold true for the other two coarse problems (fT^ and fl2^ . 



3.5. Error estimates in terms of the local energy captured. The following theo- 
rem gives and error estimates with respect to the number of eigenvectors used. Therefore, 
we obtain convergence to the reference solution when we add more and more eigenvec- 
tors to the coarse space. The error estimates are written in terms of the amount of local 
energy (of the reference solution) that is captured using the selected eigenmodes. 

Theorem 3.5. Let u]^ and u^^'^^ he the solution of Problem ^ and / flgj) . respectively. 
Put u\ = {u*}fLi. In each suhdomain we can write 



,*\2 



and ai{u*,u*) = ^Xiyt{u*f 
e=i 1=1 

where c^{u*) = mi{u*^ipl{) + mf{u*,ip-^^). The following error estimate holds 



N N, 



If* msJI * ms,II\ ^ r~i \ T t ^ I \ ^ \ ^ \II / *\2 

-Uh' .U^-Uh' ) < Ci 1 + ■ -jj 2^ 2^ \iCl[Ui) 



Proof. Using the truncated expansion of solutions, define the interpolation I^{u^) by 

= {/f «)}f=i where /f «) = f^ce{uMi 



where Ci{u*) = mi{u*, ipl^f) + mf{u*, ip-^^), i = 1, . . . ,N. Note that lf{u*) is the projec- 
tion of u* into the space spanned by the first Lj eigenvectors. Now we take v = I^{uD 
in Lemma 13.41 to obtain 

(27) dh{ul - Uh''^^, K - Uh"'") < Cidh{ul - l"{ul),ul - l"{ul)) = dh{e, e) 
where we have defined e = {ci}^^ where 

Ni 



e=Li+i 

Now we bound dh{e, e). First we observe that 

N 

(28) 4(e,e) = ^ {ai{ei,ei) +Pi{e,e)). 



i=l 
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The second term in this sum can be bounded as follows. We write 
Pi{e,e) = ^ -^-^ f Kij{ej - Cifds 



< 2 ^ — — / Kij{ejfds + 2ml{ei,ei). 



^ij(Ci) ds 



N ^ 1 6 f ^ 

i=i i=i E^cani '^'^ -'^-3 i=i 



Adding over all sub domains we get 

N N ^ ^ ^ N 

N 

(29) = 4^m^(e„e,). 

1=1 

On the other hand, if we use the increasing order of eigenvalues of eigenvalue problem 
(HZD, we get, 

(30) m\{ei, Ci) < mi{ei, e^) + (ci, e^) < -j^ — ai{ei, e^) 
which, together with (!28l) and (!29l) . gives 

N N 

Y {ai{ei, ei) + pi{e, e)) < ^(ai(ei, e^) + 4mf (e^, e^)) 

1=1 i=l 

i=i \ "^iM+iJ e=L,+i 

This completes the proof. □ 

Remark 3.6. Using our analysis, in order to obtain further hounds for the error we have 
to study the convergence of the sum \l\cii{ulY (that is, the decay of the coefficients 
ce{u*Y with increasing i). This can depend on the smoothness of the solution and it will 
be matter of further research. 

Remark 3.7. A similar result holds for the method constructed with snapshot space 
presented in Subsection \3.3[ In this case the reference solution is the solution obtaining 
by a Galerkin projection on the snapshot space. See Remark \3.1\ 
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4. Numerical experiments 

In this section we present representative numerical experiments. In particular, we 
compute the coarse (or upscaled) solution and study the error with respect to the ref- 
erence solution (or the fine-grid solution of ([3])). We choose 6 = 4 for all the numerical 
test presented here. We note that the solution of ([3]) depends on both fine-scale and 
coarse-scale parameters, h and H. We are interested mainly on the convergence (to the 
reference solution) when we sequentially add more and more basis functions. We study 
the error behavior due to the addition of coarse basis functions for fixed value of h and 
H. 

We consider the domain Q = (0, 1)^ and divide Q into N = M x M square coarse 
blocks, {fli}fLi, which are unions of fine elements. In this case H = 1/M is the coarse 
mesh parameter. Inside each subdomain Qi we generate a structured triangulation with 
m subintervals in each coordinate direction (and thus h = l/(Mm) is the fine mesh 
parameter). We consider the solution of Equation ([3]) with / = 1 and a high contrast 
coefficient described in Figuredl This coefficient is one in the white background and value 
1] in the gray regions representing high-contrast channels and high-contrast inclusions. 
Thus, 7] represents the contrast of the media, namely the ratio of the maximum and 
minimum values of k,{x). 




Figure 1. High contrast coefficient. 

In the following we compute the norm of the error e = u*^ — Uh between the fine-scale 
solution obtained by solving ([3]) and the coarse-grid solution uh, which is one of the 
following coarse-grid function: 1. u^^'^ solution of f|T6|) . 2. u^"^'^^ the solution of f lT9|) . or 
3. u^^' the solution of fl23|) . The total error is ||e||^ where || ■ ||^ ^ defined in f lTT]) . The 
relative error is computed as ||e||^ The error is divided into two quantities: 

• Interior Error: (square of the) broken if ^— semi- norm of the error 

N N ^ 

^ai(e,e) = ^ || /t/Vci 11^2(0,) . 

i=l i=l 
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Interface Error: (square of the) L^— norm of the jump of the error across the 
edges 

TV 



\2 



• Energy error: (square of the) DG bihnear form, that is, (e,e). 

4.1. Fine-grid snapshot space and original eigenvalue problem. In this Subsec- 
tion we present the numerical experiments for the method introduced in Subsection 13.11 
and show the error obtained when the dimension of the coarse space is increased. 

First, we recall that for high-contrast problems we include the eigenvectors corre- 
sponding to small eigenvalues (that asymptotically vanish as the contrast increases). 
We denote by Lf^""^^ the number of these small eigenvalues in To see the effect of 
adding more basis functions, we select additional L"'^'^ eigenvalues so the total number 
of eigenvalues selected in the block Vti is Lf^""^^ + L^'''^. We show that the error decays 
as Ll'^ increases. For the coefficient k{x) and coarse mesh shown on Figure [1] there 
is only one such and eigenvalue in each coarse-grid block VLi and therefore Lf^°'^^ = 1, 
i = 1,2, . . . , N. Figure [2] illustrates the effect of using increasing number of eigenvectors 
in the solution. We show the fine-scale solution and coarse-scale solutions computed 
with three different coarse spaces L"*^*^ = 0, 2, 11. 

Table 1. Numerical results for the Method I (see Section [3.ip . Here, 
h = 1/100, H = 1/10, fi = 10^ and n = 10^ (in parenthesis) 

Dim. Interface error Interior error Total error Amin 

100 (100) 0.026 (0.026) 0.326 (0.326) 0.3522 (0.3522) 689.4 ( 689.3) 

2 300 (300) 0.031 (0.032) 0.221 (0.220) 0.2523 (0.2516) 1562.2 (1561.7) 

4 500 (500) 0.028 (0.029) 0.171 (0.170) 0.1991 (0.1984) 2607.5 (2607.0) 

6 700 (700) 0.027 (0.027) 0.133 (0.131) 0.1600 (0.1581) 5199.4 (5199.3) 

8 900 (900) 0.026 (0.027) 0.114 (0.113) 0.1399 (0.1392) 7237.9 (7237.6) 

10 1100 (1100) 0.025 (0.025) 0.104 (0.103) 0.1293 (0.1283) 9509.1 (9509.0) 



The results for the computation of interior and interface errors are presented in Table 
[T]for h = 1/100 and H = 1/10 and two different contrasts. The convergence with respect 
to the minimum left out eigenvalue is shown in Figured (left). In this case, we solve a 
100 X 100 eigenvalue problem in each coarse block. From the results we see convergence 
to the reference solution (fine-grid solution). We also observe error decay proportional 
to the minimum left out eigenvalue across all coarse blocks. In particular, for if = 1/10, 
we need only 3 or 4 additional functions to get an interior error or the order of 17%. 
This error is computed with respect to the fine-grid solution with fine-grid parameter 
h = 1/100. Note that, in this case, we have total of four or five basis functions per 
subdomain which is comparable to the number of degrees of freedom of a classical DG 
method on the coarse grid. 
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Figure 2. Illustration of fine and coarse solutions. Fine-scale solution 
ul (top left). Coarse-scale solution u^'^'^ with L'^'^'^ = eigenvalues (top 
right). Coarse-scale solution m^*' with Lf'^'^ = 3 eigenvalues (bottom left). 
Coarse-scale solution m^*'^ with L?*^*^ = 11 eigenvalues (bottom right). 









^"Contrast 10'' 
'-'Contrast 10* 







































■ 0.5 1 1.5 X 10"' 

Figure 3. Total error of Method I (in Section [3.11) vs. 1/Amj„. Here 
h = 1/100 and the contrast is 10^. Coarse mesh size H = 1/10 and 
contrast 10^ and 10^. 

4.2. Error vs. coarse problem penalty scaling. Now we test the error when we 
change the scaling of the penalty. 
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We recall that the fine-scale problem in (Mi) uses a penalty term scaled hj 6-r--r-. In 
the classical SIPG formulation on the coarse grid one uses a penalty scaled by Sj-jj. 
Here we experiment by computing the coarse solutions with several penalties in the 
range from to 1//;, to identify a good penalty parameter for the coarse problem. 




Figure 4. The error for values of L'^ = 0, 1, . . . , 10 and the two different 
penalty scalings: | and jj. Here h = 1/100, H = 1/10 and rj = W^. 




Figure 5. The error for Method I (in Section ISTTl vs penalty scaling for 
different values of Lf'^ . Here h = 1/100, H = 1/10 and t] = 10^. The 
interior error (left) and the boundary error (right). 

For this experiment we set the contrast rj = 10^, M = 10, rii = 10, i = 1,...,N 
(and thus H = 1/10 and h = 1/100). Then, recalling that 6 = 4 for the numerical 
experiments, we have = 40 and Sj^ = 400. For these two choices of the penalty 
in Figure |4] we show the decay of the interior and interface error when adding more 
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eigenfunctions. We observe a reduction of the error as we use more and more additional 
coarse-grid basis functions. Also, we observe that the interior error are of comparable 
size (with either scaling) and that the interface error is bigger if we use the coarse penalty 
scaling. 

We fix the number of additional eigenvectors L"'*'^ and compute the interior and bound- 
ary errors when the coarse solution is computed with different penalties. We note that 
the fine-scale solution is computed with the fine-grid scaling of the penalty bilinear 
form. We repeat this experiments with Lf = 0, 1, . . . , 10. Results are shown on Figure 
13 From these results we observe that the boundary error for the coarse-grid solution 
is more sensitive to the variations in the scaling of the penalty. Indeed, for example 
M = 10, with L'^'^'^ = 4, i = 1, . . . , N, the optimal penalty coefficient is approximately 
70. These set of experiments demonstrates that one needs to properly choose the penalty 
parameter in order to balance boundary and interior errors. 

4.3. Fine-grid snapshot space and amended eigenvalue problem. Here we con- 
sider the method introduced in Subsection 13.21 We repeat the experiment described 
in Subsection 14.11 The results are displayed in Table |2] for contrast rj = 10^. Similar 

Table 2. Numerical results for Method II (Subsection 13. 2p with increas- 
ing dimension of the coarse space, h = 1/100, H = 1/10, and r] = 10^. 



j-^add 


Dim 


Interface 


Interior 


Total 


Energy 


•^min 





100 


0.0318 


0.2854 


0.3172 


0.3172 


0.0528 


2 


300 


0.0342 


0.1679 


0.2020 


0.2646 


0.0933 


4 


500 


0.0257 


0.1066 


0.1323 


0.1793 


0.1165 


6 


700 


0.0214 


0.0800 


0.1014 


0.1444 


0.2459 


8 


900 


0.0194 


0.0581 


0.0775 


0.1164 


0.3514 


10 


1100 


0.0185 


0.0566 


0.0751 


0.1132 


0.4551 



results where observed for higher contrast. We observe error reduction when next eigen- 
functions are added. Note that the results obtained by using the amended eigenvalue 
problem in Subsection 13.21 are slightly better. In this case our numerical results verify 
our theoretical error estimates in Theorem 13. 5[ Note that we report the energy error. 

4.4. Local solutions as snapshot space. Next, we consider the snapshot space given 
as in Subsection 13.31 that consists of aj— harmonic functions defined in fl2U]) . We note 
that this space of snapshots are used in the generalized multiscale finite element method 
for wave equation in [TU]. The objective of presenting these results is to show that our 
proposed DG method is flexible and one can use various snapshot spaces. 

In Table |3l we present the numerical results with contrast rj = 10^. Note, that we 
have obtained similar results for contrast t] = 10^ (not reported here). In this example, 
we choose if = 1/10 and the same setup as in the previous section. From these results 
we observe convergence to the fine-grid solution of ([3]) . We also observe that the error is 
inversely proportional to the A^m (the minimum left out eigenvalue). Note that, we are 
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Table 3. Numerical results for Method III (see Subsection 13.31) with in- 
creasing dimension of the coarse space, h = 1/100, H = 1/10 and r] = 10^. 





Dim 


Interface 


Interior 


Total 


Energy 







100 


0.032 


0.285 


0.317 


0.373 


0.0528 


2 


300 


0.034 


0.168 


0.202 


0.265 


0.0933 


4 


500 


0.026 


0.107 


0.133 


0.181 


0.1165 


6 


700 


0.021 


0.079 


0.101 


0.144 


0.2459 


8 


900 


0.019 


0.058 


0.077 


0.116 


0.3514 


10 


1100 


0.018 


0.057 


0.075 


0.114 


0.4552 



reporting the error not with respect to the reference solution (on the snapshot space) 
but with respect to the fine-grid solution. Nevertheless, we observe good results. 

We recall that we need to use the bilinear form in the eigenvalue problem in order 
to obtain error estimates. We also observe convergence in the numerical tests if we use 
rrii bilinear form (instead of mf) in the snapshot space of harmonic functions. These 
results are reported in Table HI 

Table 4. Numerical results for snapshot space and bilinear form rrii (in- 
stead of mf). Here h = 1/100, H = 1/10, and r] = 10^. 





Dim. 


Interface 


Interior 


Total 


2 


300 


0.0301 


0.1750 


0.2051 


4 


500 


0.0278 


0.1104 


0.1382 


6 


700 


0.0257 


0.0905 


0.1162 


8 


900 


0.0243 


0.0763 


0.1007 


10 


1100 


0.0223 


0.0663 


0.0986 



5. Discussions on the convergence 

In this section, we discuss the convergence of the proposed discontinuous multiscale 
finite element method. For the multiscale method in Section I3.2[ we derived error 
estimates in Theorem 13.51 Similar result holds for the coarse space described in Section 
13.31 Despite of the fact we do not write an error estimate for the multiscale space in 
Section 13.1^ which uses simple weighted eigenvalue problems, we verified convergence in 
the numerical experiments. Moreover, we observe that the interface error due to penalty 
term is smaller than the interior error calculated with the energy norm (see Table E]). 
This indicates that once we have sufficient number of multiscale basis functions per 
coarse region (that are selected properly with our spectral problem) the interface error 
is dominated by the interior error. Consequently, one can argue that the local basis 
function construction should mostly take into account the approximation property within 
coarse regions the idea that we follow in this paper. The proposed eigenvalue problem 
attempts eliminating some interior degrees of freedom via an (interior and a boundary) 
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mass matrix as discussed in our previous papers, e.g. [151 120]. By selecting mass matrix 
carefully in the local spectral problem, we represent piecewise constant functions within 
high-conductivity inclusions as the exact solution becomes almost constant within these 
regions. Note that without the k (or Kjj) weight in the mass bilinear form, we will 
be selecting large number of important modes (all fine-grid degrees of freedom within 
high- conductivity regions according to previous studies). 

We have amended the local spectral problem with a mass term on the boundary (see 
Section 13.21 and 13. 3p in order to obtain error estimates in terms of the local energy 
captured by the local space. A number of other similar eigenvalue problems that we 
tried, including adding some of the penalty terms in the eigenvalue problem, produced 
less satisfactory results. 

As we observe from our numerical results that the penalty does not need to be very 
large. Even though the boundary error decreases, it is much lower than the interior 
error for relatively small penalty terms. Thus, if a small penalty term is preferred (e.g., 
for time-dependent problems), then one can attempt to find an optimal penalty. For 
the steady state problem considered in this paper, we do not investigate this problem. 
We believe that other appropriate techniques for coupling discontinuous spectral basis 
functions should be also considered. For instance, formulations that avoid the use of 
artificial penalty terms and use approaches such as hybridized discontinuous Galerkin 
methods, [TT], or mortar methods, [B]. This is object of current research. 

One can choose other multiscale spaces such that to achieve higher accuracy with lower 
degrees of freedom. For example, as we show that using proper eigenvalue problem one 
can improve the error. Error to a reference solution can also be improved by selecting 
appropriate snapshot spaces. In general, the choice of snapshot space depends on the 
input space as discussed in Introduction (see also [17] and the local solution space which 
this input space gives). Here, we do not explore the choice of snapshot spaces since 
our goal is to show that Discontinuous Galerkin provides a nice framework to couple 
discontinuous basis functions computed locally. 
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